function emr_export(emrfiles,exportfile,endian,blockfilemodifier,cheborder,slicetiming)
% emr_export
%
%    emr_export(emrfiles,exportfile)
%    emr_export(emrfiles,exportfile,endian)
%    emr_export(emrfiles,exportfile,endian,blockfilemodifier)
%    emr_export(emrfiles,exportfile,endian,blockfilemodifier,cheborder)
%
% inputs:
% 
%           emrfiles : 
%         exportfile :
%             endian : 'b' or 'l' or 'n'
%  blockfilemodifier : If 'multiblock' is set in the EMR,
%                      then this var replaces %? (e.g. %i) in emr.efile and emr.mrfile
%          cheborder : order of the cheby projection. 0 to disable (default=0).
%        slicetiming : optional. 1 to enable slice timing (default).
%
% Author: Mads Dyrholm

if nargin<2
  exportfile = 'D.mad'
end

if nargin<3
  endian = 'n'
end

if nargin<5, cheborder = 0; end
if nargin<6, slicetiming = 1; end

% get data
emr = emr_read(emrfiles);
if isfield(emr,'multiblock')
  if (isfield(emr,'efile')), emr.efile  = sprintf(emr.efile,blockfilemodifier); end
  if (isfield(emr,'eventfile')), emr.eventfile  = sprintf(emr.eventfile,blockfilemodifier); end
  emr.mrfile = sprintf(emr.mrfile,blockfilemodifier);
end
fprintf('Importing events\n');
[epotimes,y] = emr_get_epotimes(emr);
fprintf('Reading data\n');
[X,slicetimes,voxelsize] = emr_get_data_and_slicetimes(emr,cheborder,slicetiming);

% open export file
[fid,msg] = fopen(exportfile,'w',endian);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write the volume header and data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% build header
%
%  int header
% 
%       1: size(volume,1)
%       2: size(volume,2)
%       3: size(volume,3)
%       4: size(volume,4)
%       5: Number of events
%
%   dbl header
%
%       1: TR in seconds
%     2-4: Voxelsize in millimeter

hdr = zeros(1,20);
for hdrelm = 1:4
  hdr(hdrelm) = size(X,hdrelm);
end
hdr(5) = length(epotimes);
fprintf('%i events.\n',hdr(5));

fprintf('Writing header, int\n');
fwrite(fid,length(hdr)*4,'int32');
fwrite(fid,hdr,'int32');
fwrite(fid,length(hdr)*4,'int32');


hdr(:) = 0;
hdr(1) = emr.TR;

hdr([2,3,4]) = voxelsize;
fprintf('Writing header, dbl\n');

fwrite(fid,length(hdr)*8,'int32');
fwrite(fid,hdr,'real*8');
fwrite(fid,length(hdr)*8,'int32');

% write the data
fprintf('Writing data\n');
elms = prod(size(X));
fwrite(fid,elms*8,'int32');
fwrite(fid,X,'real*8');
fwrite(fid,elms*8,'int32');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write the slice times
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Writing slice times\n');
fwrite(fid,length(slicetimes)*8,'int32');
fwrite(fid,slicetimes,'real*8');
fwrite(fid,length(slicetimes)*8,'int32');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write the events
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fprintf('Writing event info: timing');
fwrite(fid,length(epotimes)*8,'int32');
fwrite(fid,epotimes,'real*8');
fwrite(fid,length(epotimes)*8,'int32');
fprintf(', labels.\n');
fwrite(fid,length(y)*8,'int32');
fwrite(fid,y,'real*8');
fwrite(fid,length(y)*8,'int32');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% close file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fclose(fid);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [epotimes,y] = emr_get_epotimes(emr)
% 

if isfield(emr,'eventfile')
    e = load([emr.epath emr.eventfile]);
    y = e(:,1);
    epotimes = e(:,2);
else
    
    % load e channel
    e = load([emr.epath emr.efile]);
    if size(e,1)>size(e,2), e = e'; end
    e = e(emr.channel,:);

    % quantize
    e = repmat(emr.quantizer',[1 length(e)]) - repmat(e,[length(emr.quantizer) 1]);
    [dum,e] = min(abs(e));

    clear dum % CLEAR

    if ~isempty(emr.Fs)
        % edge design
        diffe = diff(e);
        f = find(diffe>0);
        try % use emr.steadyfor
            f_all = f;
            f = [];
            for i_idx = 1:length(f_all)
                if (diffe(f_all(i_idx)+(1:emr.steadyfor))==0)
                    f = [f f_all(i_idx)+1]; % +1 because of diff
                end
            end
        end
        %f = find(diff(e)>0);
        ic = ismember(e(f+1),emr.contrast);
        iv = ismember(e(f+1),emr.versus);

        clear e % CLEAR

        tc = (f(ic)+1)/emr.Fs;
        tv = (f(iv)+1)/emr.Fs;
    else
        % ... block design if Fs = []
        ic = ismember(e,emr.contrast);
        iv = ismember(e,emr.versus);

        clear e % CLEAR

        tc = (find(ic)-1) * emr.TR;
        tv = (find(iv)-1) * emr.TR;
    end

    [epotimes,cv] = sort([tc tv]);
    y = double(cv<=length(tc));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [data,slicetimes,voxelsize] = emr_get_data_and_slicetimes(emr,cheborder,slicetiming)
% 

% load datafile
nii = load_nii([emr.mrpath emr.mrfile]);
voxelsize = nii.hdr.dime.pixdim(2:4);
dim1 = size(nii.img);
dim1 = dim1(1:3);

% slice timing
if isfield(emr,'trgap') trgap = emr.trgap; else trgap = 0; end
slicetimes0 = linspace(0,emr.TR-trgap,dim1(3)+1); slicetimes0 = slicetimes0(1:end-1);
slicetimes = slicetimes0;
slicetimes(1:2:end) = slicetimes0(1:length(slicetimes0(1:2:end)));
slicetimes(2:2:end) = slicetimes0(length(slicetimes0(1:2:end))+1:end);
%%%%%

if slicetiming==0
   slicetimes = emr.TR/2 + 0*slicetimes;
end

% ROI
if isfield(emr,'roicorner1') % old version, 1 roi allowed
  data = nii.img(min(emr.roicorner1(1),emr.roicorner2(1)):max(emr.roicorner1(1),emr.roicorner2(1)),...
		 min(emr.roicorner1(2),emr.roicorner2(2)):max(emr.roicorner1(2),emr.roicorner2(2)),...
		 min(emr.roicorner1(3),emr.roicorner2(3)):max(emr.roicorner1(3),emr.roicorner2(3)),:);
  slicetimes = slicetimes(min(emr.roicorner1(3),emr.roicorner2(3)):max(emr.roicorner1(3),emr.roicorner2(3)));
elseif isfield(emr,'roisize') % new version, multiple roi's allowed, but same size
  data = [];
  slicetimes_save = slicetimes;
  slicetimes = [];
  for roinum = 1:length(emr.roicorner)
    xidx = emr.roicorner{roinum}(1)-1+(1:emr.roisize(1));
    yidx = emr.roicorner{roinum}(2)-1+(1:emr.roisize(2));
    zidx = emr.roicorner{roinum}(3)-1+(1:emr.roisize(3));
    data = cat(3,data,nii.img(xidx,yidx,zidx,:));
    slicetimes = cat(1,slicetimes,slicetimes_save(zidx));
  end
else
  data = nii.img;
end

if cheborder > 0
  dimch = size(data); %nii.hdr.dime.dim(2:5);
  % cheby boy
  %cheborder = 4;
  kind = 1;
  C = [];
  for n=0:cheborder
    C = [C;ortho_poly(1,linspace(-1,1,dimch(4)),n)];
  end

  % project
  X = double(reshape(data,[prod(dimch(1:3)) , dimch(4)]));
  X = X-X*pinv(C)*C;
  X = reshape(X,dimch);
  data = X;
end



function pl=ortho_poly(kf,x,n)

% function pl=ortho_poly(kf,x,n)
%       ==========================================================
%       Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
%                or Ln(x) or Hn(x), and their derivatives
%       Input :  KF --- Function code
%                       KF=1 for Chebyshev polynomial (First kind) Tn(x)
%                       KF=2 for Chebyshev polynomial (Second kind) Un(x)
%                n ---  Order of orthogonal polynomials
%                x ---  Argument of orthogonal polynomials
%       Output:  PL(n) --- Tn(x) or Un(x) or Ln(x) or Hn(x)
%                DPL(n)--- Tn'(x) or Un'(x) or Ln'(x) or Hn'(x)
%       =========================================================

% The only improvement in this program is it accepts vector arguments for x

% make sure that x is a row or column vector and not a matrix.
[r,c]=size(x);
if r==1 | c==1
    rowvec = 0;
    if r==1
        x=x';
        rowvec = 1;
    end
else
    error('x must be a vector, and cannot be a matrix');
end
lenx = length(x);

if n==0
    if rowvec
        pl = ones(1,lenx);
    else
       pl = ones(lenx,1);
    end
else
    pl = zeros(lenx,n);
    
    a=2;
    b=0;
    c=1;
    y0=1;
    y1=2.*x;
    
    % the i'th position in pl corresponds to the i'th term
    % don't bother storing pl = 1;
    
    pl(:,1)=2.*x;
    
    if (kf == 1)
        y1=x;
        pl(:,1)=y1;
    end
    
    for  k=2:n
        yn=(a.*x+b).*y1-c*y0;
        pl(:,k)=yn;
        y0=y1;
        y1=yn;
    end
    if rowvec
        pl = pl(:,n)';
    else
       pl = pl(:,n);
    end
end
